Analysis of a fully packed loop model arising in a magnetic Coulomb phase 
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The Coulomb phase of spin ice, and indeed the Jc phase of water ice, naturally reahse a fully- 
packed two-colour loop model in three dimensions. We present a detailed analysis of the statistics of 
these loops, which avoid themselves and other loops of the same colour, and contrast their behaviour 
to an analogous two-dimensional model. The properties of another extended degree of freedom are 
also addressed, flux lines of the emergent gauge field of the Coulomb phase, which appear as "Dirac 
strings" in spin ice. We mention implications of these results for related models, and experiments. 



Introduction: In the study of magnetism, we are nat- 
urally led to consider local degrees of freedom and their 
correlations, such as the long-range order of the local spin 
orientation in a ferromagnet. In other settings, the fun- 
damental degrees of freedom are extended, with polymers 
presenting perhaps the most familiar instance. 

In magnets in two dimensions (2D), linearly extended 
objects do occur in the form of magnetic domain walls. 
The study of such 'non-local' degrees of freedom, involv- 
ing questions such as: "What is the probability for two 
bonds to be on the same domain wall?", is related to 
problems like the geometry of the hull of a percolating 
cluster. Some beautiful theories have been developed in 
this context [H-Q. The same questions in higher dimen- 
sion do not lend themselves to similarly exact approaches, 
but fractal extended objects in 3D have been studied in 
the contexts of polymer physics Q , cosmic strings [1, @] , 
magnetic filaments in manganite materials 0, and laser 
speckles @. 

Here, we discuss a 3D frustrated magnetic system, spin 
ice [^, exhibiting two distinct extended degrees of free- 
dom, which we call loops and worms. Spin ice is un- 
usual in that its low-temperature magnetic state is nei- 
ther ordered (as in a ferromagnet) nor disordered (as in 
a paramagnet). Rather, this state is a Coulomb phase, 
where an emergent conservation law leads to algebraic 
spin correlations at low temperatures [lol |. This has re- 
cently been confirmed experimentally for the compound 
Ho2Ti207 The loops define a two-color fully packed 



loop model on the diamond lattice (the "premedial" lat- 
tice of the pyrochlore [12|), and they also appear in mod- 
els for pyrochlore compounds with two species of mag- 
netic ions I12h14|| or itinerant electrons subject to dou- 
ble exchange [15|. Worms, named after their appearance 
in Monte Carlo worm algorithms [l6j . only come in one 
fiavour and can pass through themselves and each other. 
They play a conceptually important role in the physics 
of spin ice as they are connected to "Dirac strings" and 
deconfinement in the Coulomb phase 17, l^]- Thus, be- 
sides their interest in the statistical mechanics of lattice 
models, we study loops and worms to elucidate the prop- 
erties of this new magnetic phase. 

In this paper, we numerically evaluate fundamental 
characteristics such as the probability distribution func- 



tion (PDF) of loop length radius of gyration and the 
probability for two sites separated by distance r to be on 
the same loop, C(r). We present analytical arguments to 
account for the observed regimes and their concomitant 
power laws, and contrast this behaviour to the analogous 
model in two dimensions, on the checkerboard (Fig. [1}, 
whose premedial lattice is the square lattice. Finally, we 
mention possible experimental signatures and touch on 
related models which naturally occur in the study of frus- 
trated magnets and multi-colored loop models jTH, [l^] • 

Loops and worms: In spin ice, classical Ising spins live 
on the sites of the pyrochlore, or equivalently, the bonds 
of the diamond lattice. At low temperature, an exponen- 
tially large number of degenerate ground state configu- 
rations is available, exp(A^S'p), where N is the number 
of spins and Sp = ^ log | is Pauling's ice entropy. These 
states obey the ice rules, stating that two of the bonds 
emanating from a diamond lattice site are occupied by up 
spins, the other two by down spins. They correspond to 
the allowed configurations of cubic ice /c, captured by the 
six- vertex model on the diamond lattice. The ensemble 
of these states provides the Coulomb phase exhibiting an 



emergent gauge field and algebraic correlations [10|, 112 1 . 

By coloring the links occupied by up and down spins 
differently, we obtain a fully packed (each link hosts a 
loop segment) two-color loop model. In contrast, a worm 
contains an alternating sequence of adjacent up and down 
spins (Fig. [T|). There are different possible constructions 
for a worm. We adopt an unbiased one in which a worm, 
having entered a tetrahedron, exits through one of the 
two opposite spins in that tetrahedron with equal prob- 
ability. A worm ends when it meets its initital site. 

We consider periodic systems of square or cubic geom- 
etry with L unit cells in each direction, so that there are 
N = ALP' and N = IQLi^ lattice sites for the checkerboard 
and pyrochlore cases respectively. The smallest possible 
loop is imin = 4 (6) for checkerboard (pyrochlore). We 
denote by {d) the average loop length and use subscripts 
c or 2d {p or 3d) for checkerboard (pyrochlore). 

Loop length distributions: Fig. [5] presents the PDF of 
the loop length obtained using the Monte Carlo worm 
algorithm [16|. For the checkerboard, the PDF has a 
single power-law behavior, P2d ~ l^'^" , with Tc = 2.14(1). 
For related models (sans ice rules), this quantity is known 
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FIG. 1: Left: Loops of two colors, blue (red) for up (down) 
spins, and a worm (dashed green, made of alternating up and 
down spins), on the checkerboard lattice. The way we have 
drawn the loops highlights the fact that they fuUy occupy the 
premedial square lattice. A worm can cross itself and retrace 
its path. Right: Spin Ice model on the pyrochlore lattice with 
in (blue) and out (red) spins. 
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FIG. 2: Probability distribution functions of loop lengths £ 
on the checkerboard (left) and pyrochlore (right) lattice, for 
different system sizes {L = 100, 400, 1000 in 2D and L = 
4, 14, 60 in 3D). Each axis is scaled by a power of L. Note 
the absence of one-parameter scaling in 3D at large I. Inset: 
Average loop length on the pyrochlore lattice, converging with 
system size to a finite value of 227.3 (dashed line). 



exactly to be = 2 + 1/7 (Ref. [l] and Refs. thereof). No 
such exact value is available for Tp. In fact, the situation 
for the 3D pyrochlore is dramatically different. There 
are now two different power-law regions: a short loop 
region where the PDF scales as ^ L'^£~^-^^'-^\ and then 
a crossover at £i ~ to a large-^ region with scaling 
~ £-0.98(3) _ second regime is due to the influence of 
winding loops, which close only after crossing the periodic 
boundaries. We have checked that the PDF of only non- 
winding loops has a single power law, - L3^-5/2. In 
2D, the loop distribution is dominated by non-winding 
loops at all £. This is reminiscent of Polya's theorem, 
that a random walk in 2D is recurrent (always returns 
to the initial point), while in 3D it is transient (a finite 
probability never to return) [l^. In both 2D and 3D, 
P{i) briefly increases just before the cutoff at large £. 

Winding and non-winding loop fractions: In 2D, the 
average number of winding loops is 1.86(1) for large 
L. In contrast, in the pyrochlore this number scales as 
^ InL. In both cases the number of non- winding 
loop scales as the total number of sites. In the large-L 
limit, the percentage of pyrochlore sites belonging to non- 
winding and winding loops are 6.310(5)% and 93.690(5)% 
respectively. Thus a finite but small portion of the system 
is covered by an extensive number of short loops, while 
most of the lattice sites belong to a few large loops. 

Probability to be on the same loop: In 2D, wc find 
C{r)^r~^/'^ , in analogy with 0]. By contrast, in 3D, the 
leading term in C{r) is a constant (Fig. i.e. sites far 
from each other have a nonzero probability to be on the 
same loop. This is another qualitatively new feature of 
the 3D case arising from the transient nature of 3D loops. 

Radius of gyration: This is defined as = 
J Yfi=i ~ where is the position of the ith site 
of the loop of length £ and (r) = X^i''*/^- This follows 



a power law R ^ £'^. We find ly = 0.573(5) in 2D, con- 
sistent with 4/7 [il. In 3D, wc find v = 0.500(5). This 
is the value for the random walk rather than the self- 
avoiding (:/sa3/5) walk, even though our loops are self- 
avoiding. This is similar in spirit to the observation for 
dense polymer solutions [4(| , that the need to avoid other 
loops counteracts the effects of self-avoidance. 
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FIG. 3: Probability for two spins to be on the same loop as a 
function of the distance between them, on pyrochlore lattices 
with L = 4, 8, 18, asymptoting to Coo = 0.2926(10) (dashed 
lines), found by an extrapolation to large system sizes (inset). 



Ideal- chain analogy: We next present an analysis 
which accounts for the above results by applying known 
results from the theory of random walks (ideal chains) to 
this setting. The probability for an ideal chain starting 
at Xo, to visit position x after £ steps, is p(xo,x;£) = 
(27r£)-3/2e-(^-''°)'/2^^ in the 3D continuum Thus, 
the probability to come back to the initial point is 
p(xo;£) = {2Tr£) ^^'^ . Summing over all possible start- 
ing positions, one obtains the following loop length PDF 
for non- winding loops: 
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where the factor l/£ compensates the arbitrariness in 
defining the starting position along the loop of size £. 

To understand the different exponent for long loops, we 
have to consider winding loops. These may be regarded 
as loops that start at a point in the "original" sample 
but end at an equivalent point in a "copy" sample, aris- 
ing from periodic boundary conditions. The copies cover 
all space. The probability for the loop reaching the equiv- 
alent point in the i-th nearest-neighbor copy at distance 
r,; = HiL, is P{i;£) = (27r£)-3/2e-'-?/2^ in 3D. The num- 
ber of copies which are i-th nearest neighbors scales as 
47r(r,;/L)2 for large i. Thus the total probability for a 
winding loop starting at the origin to be of size £ is 
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after approximating the sum by an integral. We note that 
p(0;£) is independent of £ and scales as . Summing 
over all possible initial points, the total probability to 
have a loop of size £ becomes 



(3) 



independent of system size L. The crossover between the 
two behaviors occurs at the length scale £i^L'^, i.e. when 
the radius of gyration of a loop reaches the system size. 

The average loop lengths {£) in both checkerboard and 
pyrochlore cases are finite. For the checkerboard, the 
PDF exponent Tc being larger than 2 ensures a non- 
diverging {£)c. With the approximation that P2d{t) 
£~'^'' at all even £ starting from £min — 4, one estimates 
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24.9, (4) 



very close to the numerically obtained average 24.68(3). 

For the pyrochlore, the combination of non-winding 
{P^J ~ L3^-^/2) and winding loops {P^i ~ ^"^) conspire 
to produce a remarkably large but finite average loop 
length: (£)p = 227.5(5). Even though they increase {£)p 
by an order of magnitude compared to the 2D model, 
winding loops do not manage to make it divergent. 

Considering non-winding and winding loops sepa- 
rately, we find that the non-winding loop average sat- 
urates with L, to the value (£)n.w. = 14.34(4). (Eq. (g]) 
adapted to the pyrochlore gives 15.3). The winding loops 
by themselves have a diverging average (^)w ~ L'^/lnL. 

Concerning the probability for two sites to be on the 
same loop, we note that the number of spin pairs belong- 
ing to the same loop of length £ is ^£{£ — 1) x P{£)^ while 



there is a total of ^N{N — 1) pairs of spins in the system. 
The probability averaged over all pairs is thus 
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P34£)d£ ~ 0(L°) (5) 



The non-winding loop contribution to this integral van- 
ishes at large L, but the winding loop contribution leads 
to a constant term Coo = C{r oo) (Fig. [3]). Indeed, we 
note that there will be several loops in a large system, 
each of which contains a finite fraction, /i, of all sites. 
(We found fo « 0.41, /i w 0.30, /a « 0.12 . . . ). 

The 'scaling function' displayed in Fig. [5] shows that 
in 3D, there are in fact two length scales, namely when 
the radius of gyration hits the system size, h ^ , and 
when the loop explores the full volume, I2 ^ . The 
behaviour for I > h obviously is influenced by the nature 
of the boundary conditions, e.g. the loops winding many 
times can break up into several loops terminating on the 
surface when open boundaries are considered. 

For 2D, we note an intriguing similarity to the geom- 
etry of percolation hulls, which have a size distribution 
scaling with the same exponent we find above poj . This 
we can rationalise as follows. Our loops can be thought of 
as defining a specific percolation problem, with checker- 
board lattice sites painted in either colour with proba- 
bility p ~ 1/2. It has been argued that neither short- 
range nor sufficiently rapidly decaying algebraic correla- 
tions between the occupied bonds influence the percola- 
tion critical exponents [2l| . The 2D spin ice correlations 
indeed decay sufficiently fast, as r"^. An additional fea- 
ture of our loop ensemble is that loops of the same colour 
cannot cross or branch. From this it follows that each 
loop is at the same time a hull, whose length/perimeter 
should scale as L^/-* (Fig.©. 

Worms: Wc now address properties of the worms. 
These are efficiently evaluated as our Monte-Carlo al- 
gorithm is in fact based on constructing worms, as flip- 
ping all spins in a worm conserves the ice-rules. Because 
worms can retrace parts of their path, their actual length 
can be longer than the system size. We therefore look 
at the distribution Q{X) of the number of spins flipped 
during a worm update (Fig. 3]), i.e., the number of sites 
visited by the worm an odd number of times. Again, two 
regimes appear in 3D, an ^"■^/^ behavior (numerical ex- 
ponent -1.48(4)) for ^<£i, and a ^-independent region for 
larger £. These (—3/2 and 0) are random-walk exponents 
that can be derived similarly as in the loop case, Eqs. ([1]) 
and ([3]). The exponents for Q{X) are shifted from the 
loop P{£) exponents (—5/2 and —1) by 1; indeed if one 
asked the question "How long is the loop on which a ran- 
domly chosen spin sits ?" , one would also get exponents 
-3/2 and 0. 

The winding worm regime implies that 3D worms on 
average visit a finite fraction of the system size (« 15%). 
Hence, a worm algorithm where a Monte Carlo update is 
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FIG. 4: Data points: probability distribution of number of 
spins flipped by worms (L — 4, 10, 40). Solid line: probability 
distribution of full worm length {L = 40) which weighs each 
site by the number of times it has been visited. This departs 
parametrically only for the longest worms, where multiple 
visits to a site become appreciable. 



made by flipping the spins in a worm only requires a finite 
(typically 5—10) constant nmnber of worms to decorrelate 
the system in 3D, independent of system size L, while 
the different behaviour of winding worms in 2D requires 
an increasing number of worm updates for decorrelation 
with L. 

Conceptual and experimental implications: For both 
2D and 3D, we have recovered exponents originally de- 
rived for systems without ice rules. The divergence-free 
conditions of Coulomb phases seem to have little effect 
on the 'universal' loop and worm statistics we have con- 
sidered, which was not a priori to be expected, as the 
ice-rules impose algebraic spin correlations. 

In the context of spin ice, worm statistics plays an 
important conceptual role as worms describe the tension- 
free (emergent) flux loops which are characteristic of the 
Coulomb phase. Here we have shown that their statistics 
resembles that of an ideal chain. Since, in a magnetic 
field, the worms have been utilised as an experimental 
diagnostic of Coulomb phase physics [l8| , a more detailed 
study of their statistics as a function of field strength 
would be worthwhile. 

Regarding the loops, these are most naturally probed 
as a non-local correlation function. For instance, for a 
system of electrons which can hop only along a loop, the 
existence of the winding loops shows up in conductivity 
properties, as we will discuss elsewhere [isj. 

To the best of our knowledge, our loops first appeared 
in Villain's seminal paper in the context of a model of 
pyrochlore Heisenberg magnets with two species of mag- 
netic ions [l2| - [l3 |. Villain already noted the possibility 
of two (recurrent and transient) loop populations, which 
also show up for cosmic strings 0, Ig] and laser speck- 
les Hi. 

In Villain's model, each loop has a distinct colour en- 
coding the direction of its Heisenberg spin (which is con- 
tinuously variable and hence the number of colours is 
infinite), and spins are correlated only if they belong 
to the same loop. A finite value of Coo thus implies 



long-range spin order. Here we comment on two note- 
worthy features. Firstly the presence of several loops of 
size 0{L^) implies coexistence of several spatially inter- 
twined but independent populations of long-range or- 
dered spins. Secondly, a short-range interacting classical 
spin Hamiltonian of the kind envisaged by Villain allows 
gapless excitations for the Heisenberg spins. Thermal 
fluctuations out of this set of states could lead to locking 
of these populations into a more conventional collinear 
ordered structure. This is clearly a productive field for 
more detailed studies. 

If one assigns a weight to each loop reflecting the num- 
ber of flavours, n, it can take (which is not what is done in 
Villain's model, whose colours are used to distinguish the 
loops, not to weigh them), it becomes more favourable to 
have many, short loops as n grows. There is an extended 
literature on loop models, with some 3D work [23, [23j . 
from which it is known that a phase transition results. In 
our case, an n cxd state is the "hexagonal protectorate" 
loop crystal proposed in the context of magnetodistortive 
phase transitions in the spinel compound ZnCr204 [24 1. 



This has only loops of length 6 and breaks lattice trans- 
lation symmetry. 

In summary, we have presented an analysis of a set of 
extended degrees of freedom arising in an exotic phase of 
a three-dimensional magnet. This we hope will be of dual 
interest both from a statistical mechanics and a magnetic 
materials relaxation perspective. Further studies of these 
models, e.g. a quantum version thereof, are ongoing [25j . 
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